using System;
using System.Collections.Generic;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.LinearAlgebra;
using Atomic.Structures;

namespace Atomic.Thermodynamics
{
	[Serializable]
	public class RedlichKisterPolynomial
	{
		private int n, m;

		public RedlichKisterPolynomial(int order, Atom[] atoms, double[] coefficients)
			: this(order, atoms, coefficients, new double[atoms.Length])
		{
		}

		public RedlichKisterPolynomial(int order, Atom[] atoms, double[] coefficients, double[] referenceEnergies)
		{
			if (order < -1)
			{
				throw new ArgumentException();
			}

			if (coefficients == null || atoms == null || referenceEnergies == null)
			{
				throw new ArgumentNullException();
			}

			Atoms = new ReadOnlyList<Atom>(atoms);
			Coefficients = new ReadOnlyList<double>(coefficients);
			ReferenceEnergies = new ReadOnlyList<double>(referenceEnergies);

			// Number of atoms.
			n = Atoms.Count;

			// Number of fittings coefficients per pair of atoms.
			m = order + 1;

			if (Coefficients.Count != m * n * (n - 1) / 2)
			{
				throw new ArgumentException("Invalid number of coefficients.");
			}

			if (ReferenceEnergies.Count != n)
			{
				throw new ArgumentException("Invalid number of reference energies.");
			}
		}

		public static RedlichKisterPolynomial Fit(int order, IEnumerable<CompositionMapping> energies)
		{
			if (energies == null)
			{
				throw new ArgumentNullException();
			}

			return Fit(order, energies.ToArray(), null, false);
		}

		public static RedlichKisterPolynomial Fit(int order, IEnumerable<CompositionMapping> energies, IEnumerable<AtomMapping> referenceEnergies)
		{
			if (energies == null || referenceEnergies == null)
			{
				throw new ArgumentNullException();
			}

			return Fit(order, energies.ToArray(), referenceEnergies.ToArray(), false);
		}

		public static RedlichKisterPolynomial FitReferenceEnergies(int order, IEnumerable<CompositionMapping> energies)
		{
			return Fit(order, energies, ExtractReferenceEnergies(energies));
		}
		
		public static RedlichKisterPolynomial FitGroundState(int order, IEnumerable<CompositionMapping> energies)
		{
			if (energies == null)
			{
				throw new ArgumentNullException();
			}

			return Fit(order, energies.ToArray(), null, true);
		}

		public static RedlichKisterPolynomial FitGroundState(int order, IEnumerable<CompositionMapping> energies, IEnumerable<AtomMapping> referenceEnergies)
		{
			if (energies == null || referenceEnergies == null)
			{
				throw new ArgumentNullException();
			}

			return Fit(order, energies.ToArray(), referenceEnergies.ToArray(), true);
		}

		public static RedlichKisterPolynomial FitGroundStateReferenceEnergies(int order, IEnumerable<CompositionMapping> energies)
		{
			return FitGroundState(order, energies, ExtractReferenceEnergies(energies));
		}

		private static RedlichKisterPolynomial Fit(int order, CompositionMapping[] energies, AtomMapping[] referenceEnergies, bool groundState)
		{
			if (order < -1)
			{
				throw new ArgumentException();
			}

			// Select group state energies only.
			CompositionMapping[] energies0 = energies.GroupBy(cm => cm.Composition).Select(g => g.OrderBy(cm => cm.Value).First()).ToArray();

			if (groundState)
			{
				energies = energies0;
			}
			else if (!groundState && energies0.Length != energies.Length)
			{
				throw new ArgumentException("Multiple energies specified for same composition.");
			}

			// Number of fittings coefficients per pair of atoms.
			int m = order + 1;

			// Concatenate arrays of atoms and determine distinct atoms.
			Atom[] atoms = energies.SelectMany(c => c.Composition.Atoms).Distinct().ToArray();
			int n = atoms.Length;

			double[] referenceEnergies0 = new double[n];
			if (referenceEnergies != null)
			{
				// Now that reference energies are specified, we require all of them to be specified.
				for (int i = 0; i < n; i++)
				{
					AtomMapping[] mappings = referenceEnergies.Where(am => am.Atom == atoms[i]).OrderBy(am => am.Value).ToArray();

					if (mappings.Length == 0)
					{
						throw new ArgumentException("Missing reference energy.");
					}

					if (!groundState && mappings.Length > 1)
					{
						throw new ArgumentException("Multiple reference energies specified.");
					}

					referenceEnergies0[i] = mappings[0].Value;
				}
			}

			// Number of fitting points.
			int u = energies.Length;

			double[] b = new double[u];
			double[,] a = new double[u, m * n * (n - 1) / 2];
			for (int t = 0; t < u; t++)
			{
				// Reference energy.
				double y = 0.0;

				for (int i = 0, l = 0; i < n; i++, l++)
				{
					double x1 = energies[t].Composition.Fraction(atoms[i]);

					// Only one term per distinct pair of atoms.
					for (int j = i + 1; j < n; j++, l++)
					{
						double x2 = energies[t].Composition.Fraction(atoms[j]);

						// Repeat from regular solution model up to specified order.
						double c = 1.0;
						for (int k = 0; k < m; k++, l++)
						{
							a[t, l] = x1 * x2 * c;
							c *= x1 - x2;
						}
					}

					// Add the reference energy term.
					y += x1 * referenceEnergies0[i];
				}

				// Target for least-square fitting. Subtract reference energy.
				b[t] = energies[t].Value - y;
			}

			Vector coefficients;
			if (m > 0)
			{
				// Least-square fitting using pseudo-inversion.
				coefficients = SingularValueDecomposition.PseudoInverse(new Matrix(a)) * new Vector(b);
			}
			else
			{
				// Allow reference energy fitting with no coefficients.
				coefficients = Vector.Zero(0);
			}

			return new RedlichKisterPolynomial(order, atoms, coefficients.ToArray(), referenceEnergies0);
		}

		public static IEnumerable<AtomMapping> ExtractReferenceEnergies(IEnumerable<CompositionMapping> energies)
		{
			if (energies == null)
			{
				throw new ArgumentNullException();
			}

			foreach (CompositionMapping energy in energies)
			{
				foreach (Constituent constituent in energy.Composition.Constituents)
				{
					if (constituent.Count == 1.0)
					{
						yield return new AtomMapping(constituent.Atom, energy.Value);
					}
				}
			}
		}

		private void CheckAtoms(Composition composition)
		{
			foreach (Atom atom in composition.Atoms)
			{
				if (atom == Atom.Undefined)
				{
					// Allow this for testing purposes.
					continue;
				}

				if (!Atoms.Contains(atom))
				{
					throw new ArgumentException("Composition contains atom not included in the fitting.");
				}
			}
		}

		public double Energy(Composition composition)
		{
			CheckAtoms(composition);

			double y = 0.0;

			for (int i = 0, l = 0; i < n; i++, l++)
			{
				double x1 = composition.Fraction(Atoms[i]);

				// Only one term per distinct pair of atoms.
				for (int j = i + 1; j < n; j++, l++)
				{
					double x2 = composition.Fraction(Atoms[j]);

					// Repeat from regular solution model up to specified order.
					double c = 1.0;
					for (int k = 0; k < m; k++, l++)
					{
						y += x1 * x2 * c * Coefficients[l];
						c *= x1 - x2;
					}
				}

				// Add the reference energy term.
				y += x1 * ReferenceEnergies[i];
			}

			return y;
		}

		public double ReferenceEnergy(Composition composition)
		{
			CheckAtoms(composition);

			double y = 0.0;

			for (int i = 0, l = 0; i < n; i++, l++)
			{
				y += composition.Fraction(Atoms[i]) * ReferenceEnergies[i];
			}

			return y;
		}

		public double PartialChemicalPotential(Composition composition, Atom atom) // Or: ChemicalPotential
		{
			CheckAtoms(composition);

			// Compute derivative with respect to the fractional amount of the specified atom. Atoms are added without removing vacancies.

			double y_x = 0.0;

			for (int i = 0, l = 0; i < n; i++, l++)
			{
				double x1 = composition.Fraction(Atoms[i]);
				double x1_x = Atoms[i] == atom ? 1.0 : 0.0;

				// Only one term per distinct pair of atoms.
				for (int j = i + 1; j < n; j++, l++)
				{
					double x2 = composition.Fraction(Atoms[j]);
					double x2_x = Atoms[j] == atom ? 1.0 : 0.0;

					// Repeat from regular solution model up to specified order.
					double c = 1.0;
					double c_x = 0.0;
					for (int k = 0; k < m; k++, l++)
					{
						y_x += ((x1_x * x2 + x1 * x2_x) * c + x1 * x2 * c_x) * Coefficients[l];
						c_x = c_x * (x1 - x2) + c * (x1_x - x2_x);
						c *= x1 - x2;
					}
				}

				// Add the reference energy term.
				y_x += x1_x * ReferenceEnergies[i];
			}

			return y_x;
		}

		public double ChemicalPotential(Composition composition, Atom atom) // Or: TotalChemicalPotential
		{
			if (!Atoms.Contains(Atom.Vacancy))
			{
				throw new Exception("Can't compute the this with no vacancies in the fitting.");
			}

			if (atom == Atom.Vacancy)
			{
				throw new Exception("Can't compute this for a vacancy.");
			}

			CheckAtoms(composition);

			// Compute derivative with respect to the fractional amount of the specified atom. Atoms are added at the cost of removing vacancies.

			double y_x = 0.0;

			for (int i = 0, l = 0; i < n; i++, l++)
			{
				double x1 = composition.Fraction(Atoms[i]);
				double x1_x = Atoms[i] == atom ? 1.0 : (Atoms[i] == Atom.Vacancy ? -1.0 : 0.0);

				// Only one term per distinct pair of atoms.
				for (int j = i + 1; j < n; j++, l++)
				{
					double x2 = composition.Fraction(Atoms[j]);
					double x2_x = Atoms[j] == atom ? 1.0 : (Atoms[j] == Atom.Vacancy ? -1.0 : 0.0);

					// Repeat from regular solution model up to specified order.
					double c = 1.0;
					double c_x = 0.0;
					for (int k = 0; k < m; k++, l++)
					{
						y_x += ((x1_x * x2 + x1 * x2_x) * c + x1 * x2 * c_x) * Coefficients[l];
						c_x = c_x * (x1 - x2) + c * (x1_x - x2_x);
						c *= x1 - x2;
					}
				}

				// Add the reference energy term.
				y_x += x1_x * ReferenceEnergies[i];
			}

			return y_x;
		}

		public ReadOnlyList<Atom> Atoms
		{
			get;
			private set;
		}

		public ReadOnlyList<double> Coefficients
		{
			get;
			private set;
		}

		public ReadOnlyList<double> ReferenceEnergies
		{
			get;
			private set;
		}
	}

	/*[Serializable]
	public class RedlichKisterPolynomial : IPlotFunction
	{
		private double y0, y1;
		private int m;
		private Vector l;

		public RedlichKisterPolynomial(int order, double[] x, double[] y)
			: this(order, x, y, 0.0, 0.0)
		{
		}

		public RedlichKisterPolynomial(int order, double[] x, double[] y, double y0, double y1)
		{
			this.y0 = y0;
			this.y1 = y1;

			// Number of fitting points.
			int n = x.Length;

			if (y.Length != n)
			{
				throw new ArgumentException();
			}

			// Number of parameters to fit.
			m = order + 1;

			double[] b = new double[n];
			double[,] a = new double[n, m];
			for (int i = 0; i < n; i++)
			{
				// Simplified 1-dimensional problem, currently.
				double x1 = x[i];
				double x2 = 1.0 - x1;

				// Target for least-square fitting. Subtract reference energy.
				double lambda = x1;
				b[i] = y[i] - ((1.0 - lambda) * y0 + lambda * y1);

				double c = 1.0;
				for (int j = 0; j < m; j++)
				{
					a[i, j] = x1 * x2 * c;
					c *= (x1 - x2);
				}
			}

			// Least-square fitting using pseudo-inversion.
			l = SingularValueDecomposition.PseudoInverse(new Matrix(a)) * new Vector(b);
		}

		public double Value(double x)
		{
			// Simplified 1-dimensional problem, currently.
			double x_A = x;
			double x_B = 1.0 - x_A;

			// Reference energies.
			double y_A = y0;
			double y_B = y1;

			// Reference energy.
			double y = (1.0 - x_A) * y_A + (1.0 - x_B) * y_B;

			double c = 1.0;
			for (int j = 0; j < m; j++)
			{
				y += x_A * x_B * l[j] * c;
				c *= x_A - x_B;
			}

			return y;
		}

		public double Derivative(double x)
		{
			// Simplified 1-dimensional problem, currently.
			double x_A = x;
			double x_A_x = 1.0;
			double x_B = 1.0 - x_A;
			double x_B_x = -x_A_x;

			// Reference energies.
			double y_A = y0;
			double y_B = y1;

			// Reference energy.
			double y = (1.0 - x_A) * y_A + (1.0 - x_B) * y_B;
			double y_x = -(x_A_x * y_A + x_B_x * y_B);

			double c = 1.0;
			double c_x = 0.0;
			for (int j = 0; j < m; j++)
			{
				y += x_A * x_B * l[j] * c;
				y_x += x_A_x * x_B * l[j] * c + x_A * x_B_x * l[j] * c + x_A * x_B * l[j] * c_x;

				c *= x_A - x_B;
				c_x *= x_A_x - x_B_x;
			}

			return y_x;
		}
	}*/
}
